Methods and devices for labeling and/or matching

ABSTRACT

Devices, such as computer readable media, and methods, such as automated methods, for labeling and/or matching. Some of the devices and methods are particularly useful for anatomical labeling of human airway trees. Some of the devices and methods are particularly useful for matching branch-points of human airway trees from represented in two or more graphs.

CROSS-REFERENCE(S) TO RELATED APPLICATION(S)

This is a continuation of co-pending U.S. application Ser. No.11/122,974, now U.S. Pat. No. 8,155,403, filed May 5, 2005, which claimspriority to U.S. Provisional Patent Application Ser. No. 60/568,184,filed May 5, 2004, the entire contents of both of which (including theappendices) are expressly incorporated by reference.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH

The government may own rights in this invention pursuant to NationalInstitute of Health grant number HL-064368.

BACKGROUND

1. Field of the Present Methods and Devices

The present methods and devices relate generally to the fields oflabeling and matching. More particular, they relate to labelinggraphical representations of trees, such as automated anatomicallabeling of human airway trees. They also relate to matchingcorresponding points of at least two graphical representations, such asautomated matching of corresponding branch-points of at least twographical representations of a tree, such as a human airway tree.

2. Description of Related Art

Lung diseases like lung cancer, emphysema, and cystic fibrosis are asignificant cause of disability and premature death in westerncountries. In North America, for example, fatalities from lung canceroutnumber those from colon, breast, and prostate cancer combined.

Lung imaging plays a crucial role in the diagnosis, study, and treatmentof lung disorders as well as in physiological studies concerned withpulmonary functionality. Modern multidetector-row CT scanners (MDCT)provide a wealth of information. Volumetric lung images of the size ofseveral hundred MBytes are not uncommon. The manual analysis of theseimages is often time-consuming, tedious, and error prone. And with thehigh volume of scans taken, manual analysis is in many cases noteconomical.

The quantitative assessment of intrathoracic airway trees is importantfor the objective evaluation of the bronchial tree structure andfunction. Functional understanding of pulmonary anatomy as well as thenatural course of respiratory diseases like asthma, emphysema, cysticfibrosis, and many others is limited by our inability to repeatedlyevaluate the same region of the lungs time after time and performaccurate and reliable positionally corresponding measurements.

Branch-point matching and anatomical labeling are both tedious anderror-prone to perform manually. Working with human in-vivo data poseschallenges. In-vivo trees deviate from ideal trees because of anatomicalvariations and because of false-branches introduced by imperfections inthe preceding segmentation and skeletonization processes.

Few attempts at automating the branch-point matching process have beenmade. Pisupati et al. (1996a) and Pisupati et al. (1996b) presented amatching algorithm based on dynamic programming, which was only appliedto very similar pairs of canine trees. Pisupati et al. stated that theyexpect the method to fail on human in-vivo scans. Park (2002) presenteda tree-matching method based on an association graph (Pelillo et al.(1999)), but his method was applied only to phantom data and does nottolerate false branches.

Publications about automated anatomical labeling are similarly sparse.Mori et al. (2000) presented a knowledge-based labeling algorithm. Theproposed algorithm was only applied to incomplete trees (about 30branches per tree), and the built-in knowledge base did not incorporateanatomical variations. Additionally, the algorithm is sensitive tomissing and added (false) branches. Kitaoka et al. (2002) developed abranch-point labeling algorithm that uses a mathematical phantom asreference. Labels are assigned by matching the target tree against thisphantom. The method cannot automatically handle false branches—they haveto be pruned manually in a preprocessing step.

Other disclosures concerning an earlier version of the present labelingmethods are described in Tschirren et al. (2003), Tschirren et al.(2002a) and Tschirren et al. (2002b).

SUMMARY

Certain embodiments of the present devices include a computer readablemedium having machine readable instructions for accessing a firstrepresentation of an airway tree of a human subject; accessing a secondrepresentation of the airway tree of the human subject; andautomatically pruning at least a portion of the first representation.

Other embodiments of the present devices, including embodiments ofcomputer systems, having additional and/or different features arediscussed below.

Certain embodiments of the present methods include an automated methodcomprising accessing a first representation of an airway tree of a humansubject; accessing a second representation of the airway tree of thehuman subject; and automatically pruning at least a portion of the firstrepresentation.

Other embodiments of the present methods having additional and/ordifferent features are discussed below.

BRIEF DESCRIPTION OF THE DRAWINGS

The following drawings illustrate by way of example and not limitation.The patent or application file contains at least one drawing executed incolor. Copies of this patent or patent application publication withcolor drawings will be provided by the Office upon request and paymentof the necessary fee.

FIG. 1A shows a human airway tree with anatomical labels that areassigned to branch-points and are based on segment names. This figure isbased on artwork published in Boyden (1955).

FIG. 1B shows a typical topological relationship of segments (variationsare possible between individuals).

FIG. 2 depicts an algorithm, which is based on a breadth-first search,that may be used to compute R_(v), discussed in greater detail below.

FIG. 3A shows a portion of a segmented, skeletonized airway tree beforepruning.

FIG. 3B shows that edge 1 in FIG. 3A is a false branch.

FIG. 3C shows the result of pruning the false branch shown in FIG. 3A,and replacing the vertex from which it extended with an edge.

FIGS. 4A and 4B show suitable nomenclature for rigid registration. Thedepicted anatomical names correspond to names in FIG. 1B.

FIG. 5A is a total view showing the result of using an embodiment of thepresent matching methods for intra-subject matching, based on an in-vivoTLC and FRC scan. The FRC tree is shown in green and the TLC tree isshown in purple. The matchings are represented by red lines.

FIG. 5B is a detail view of the tree pair shown in FIG. 4A. The falsebranches near the top center part of the purple tree do not negativelyaffect the matching result.

FIG. 6 shows an example of the notation that can be used in thepopulation average described below.

FIGS. 7A and 7B show variation of single segment measures in apopulation average. FIG. 7A shows the mean ∀ standard deviation ofsegment length. The trachea was not included because its length variesdepending on the start-point of the scan. FIG. 7B shows the standarddeviation σ_(φ) of angle φ measured from the average spatial orientationof the segment and the spatial orientation of the segment in theindividual trees. Note that σ_(φ)=0 for EndLMB_s because of rigidregistration (EndLMB_s is always aligned, in this embodiment, with they-axis.).

FIGS. 8A-8C show the variation of various inter-segment measurescompiled from measurements made on 17 different trees. FIG. 8A showsvariability of difference vector between any two anatomical segments,and standard deviation of the angle between the mean difference vectorand the difference vector for the individual trees. FIG. 8B showsstandard deviation of angle between direction vectors of any twoanatomical segments. Note the relatively constant and moderate values.FIG. 8C shows mean (red) ∀ standard deviation (blue) of segment lengthratio for every pair of segments (without trachea). This figureillustrates the great variability and hence the limited value of thismeasure to embodiments of the present labeling methods. Note that thered line is centered between the blue lines. The apparent asymmetry is avisual illusion.

FIG. 9 illustrates an example of a long “false” branch.

FIG. 10 illustrates one manner of addressing a suspected false branchthrough the introduction of a “parallel” edge. In this figure, edge e₃might be a false branch. An edge e_(p) parallel to edges e₁ and e₂ isadded.

FIGS. 11A-11C illustrate a sequence of labeling steps designed to reducecomputing time by splitting the task into sub-problems. Segments labeledduring a step are gray. The labels ‘s’ and ‘t’ mark start- andterminal-segments, respectively.

FIG. 12 shows an example of the results of branch-point labeling usingone embodiment of the present labeling methods.

FIG. 13 illustrates an example of a terminal-label problem. It is notentirely in this figure which of the two branch-points is the trueend-point of the RB10 segment.

DESCRIPTION OF ILLUSTRATIVE EMBODIMENTS

The terms “comprise” (and any form of comprise, such as “comprises” and“comprising”), “have” (and any form of have, such as “has” and“having”), “include” (and any form of include, such as “includes” and“including”) and “contain” (and any form of contain, such as “contains”and “containing”) are open-ended linking verbs. Thus, an automatedmethod “comprising” accessing a first representation of an airway treeof a human subject; accessing a second representation of the airway treeof the human subject; and (c) automatically pruning at least a portionof the first representation is a method that includes at least theserecited steps, but is not limited to only possessing these recitedsteps.

Similarly, a computer readable medium “comprising” machine readableinstructions for accessing a first representation of an airway tree of ahuman subject; accessing a second representation of the airway tree ofthe human subject; and automatically pruning at least a portion of thefirst representation is a computer readable medium that has machinereadable instructions for implementing at least these recited steps, butalso covers media having machine readable instructions for implementingadditional, unrecited steps.

The terms “a” and “an” are defined as one or more than one, unless thisapplication expressly requires otherwise. The term “another” is definedas at least a second or more.

Descriptions of well known processing techniques, components andequipment are omitted so as not to unnecessarily obscure the presentmethods and devices in unnecessary detail. The descriptions of thepresent methods and devices, including those in the appendices, areexemplary and non-limiting. Certain substitutions, modifications,additions and/or rearrangements falling within the scope of the claims,but not explicitly listed in this disclosure, may become apparent tothose or ordinary skill in the art based on this disclosure.

Human airway trees show a certain degree of similarity when comparedacross multiple subjects. The similarities can be observed up to theentry-points into the sub-lobes. This includes all 33 anatomically namedsegments (see FIG. 1A) that are commonly used in pulmonary physiologyBoyden (1955). Beyond that, the branching pattern becomes more randomand every human exhibits a different branching pattern, much likefingerprints differ between individuals.

The number of airway-tree branch-points that can be matched depends onthe source of the two trees to be compared. In inter-subject matching(trees originating from different subjects), only branch-pointsassociated with the anatomically named segments can be matched. Inintra-matching (trees originating from the same subject, e.g., imaged atdifferent lung volumes, or at different points of time), it is possibleto match points beyond the anatomically named ones.

Prior to either matching or labeling, an airway tree may be segmentedusing any suitable technique from suitable data. For example,segmentation may be accomplished on volumetric computed tomography (CT)images using processes such as those described in, e.g., (a) Chapter 2of Appendix 1 of the provisional application to which this applicationclaims priority and incorporates by reference, or (b) Appendix 4 of theprovisional application to which this application claims priority andincorporates by reference (to the extent the two descriptions differ).The segmented airway tree may then be skeletonized using any suitabletechnique, such as the technique as disclosed, for example, in Palágyiet al. (2003).

Branch-points may then be detected in the skeletonization result and thetree may be represented as a directed acyclic graph (DAG). This DAG maybe used as input data for both present branch-point matching andanatomical labeling methods.

In one respect, the present devices may be characterized as a computerreadable medium comprising machine readable instructions for: (a)accessing a first representation of an airway tree of a human subject;(b) accessing a second representation of the airway tree of the humansubject; and (c) automatically pruning at least a portion of the firstrepresentation. In one respect, the present methods may be characterizedas an automated method comprising: (a) accessing a first representationof an airway tree of a human subject; (b) accessing a secondrepresentation of the airway tree of the human subject; and (c)automatically pruning at least a portion of the first representation.The representations that are accessed consistent with such devices andmethods may take the form of data.

In another respect, the present devices may be characterized as acomputer readable medium comprising machine readable instructions for:(a) accessing a first representation, the first representation includinga segmented, skeletonized representation of an airway tree at a firstlung volume of a human subject; (b) accessing a second representation,the second representation including a segmented, skeletonizedrepresentation of the airway tree at a second lung volume of a humansubject; and (c) automatically pruning at least a portion of the firstrepresentation. In another respect, the present methods may becharacterized as an automated method comprising: (a) accessing a firstrepresentation, the first representation including a segmented,skeletonized representation of an airway tree at a first lung volume ofa human subject; (b) accessing a second representation, the secondrepresentation including a segmented, skeletonized representation of theairway tree at a second lung volume of a human subject; and (c)automatically pruning at least a portion of the first representation.

In one respect, embodiments of the present branch-point matching methodsmay be characterized as having the following steps, which can beperformed automatically: delete (prune) spurious branches from inputtrees; align the input trees by performing a rigid registration; findand match major branch-points; and match sub-trees underneath majorbranch-points, one pair of sub-trees at a time.

In another respect, embodiments of the present branch-point matchingmethods may be characterized as having the following steps, which can beperformed automatically: delete (prune) spurious branches from inputtrees; find major branch-points (e.g., the carina and the end-points ofthe two main bronchi); rigid registration based on major branch pointsfound (input trees are aligned based on the major branch-points found inthe previous step); find remaining major branch-points; match remainingmajor branch-points; and match sub-trees (remaining branch-pointcorrespondences may be found by matching, e.g., only one pair ofsub-trees at a time).

Embodiments of both the present branch-point matching and the anatomicallabeling methods may be based on matching hierarchical structures(mathematical graphs) using association graphs Pelillo (1999). Usingassociation graphs for finding graph isomorphisms is a well knowntechnique that was introduced more than 20 years ago Ballard and Brown(1982) refined the method, allowing it be used on hierarchicalstructures, such that topological and inheritance relationships werepreserved.

An association graph is an auxiliary graph structure derived from thetwo graph structures to be matched. A graph G=(V, E) comprises a set ofvertices V and a set of edges E. For two graphs G₁ and G₂, anassociation graph G_(ag)=(V_(ag), E_(ag)) consists of the verticesV_(ag)=V₁×V₂. Thus, it contains a vertex for every possible pair ofvertices in G₁ and G₂. Two vertices in G_(ag) are connected with an edgeif and only if the corresponding vertices in G₁ and G₂ stand in the samerelationship to each other (e.g., inheritance relationship, topologicaldistance, etc.). The maximum clique in the association graph (thebiggest sub-set of vertices where every pair of vertices is connected byan edge) corresponds to the maximal subtree isomorphism. Each vertexcontained in the maximum clique represents a pair of matching verticesin the original graphs. It has been shown (Garey and Johnson, 1979) thatthe problem of finding the maximum clique is an NP-complete problem. Nomethod is known to exist that can efficiently find the maximum clique ina graph with a high number of vertices and edges. However, embodimentsof the present methods may reduce the time required to match two normalsized in-vivo trees to about 1 to 3 seconds.

Pelillo et al. (1999) use path length and level difference (generationnumber) as a measure of equivalence when adding edges to the associationgraph. This is “hard” measure in the sense that only topologicallyidentical sub-trees can be matched. In contrast, embodiments of thepresent methods use a measure of equivalence that incorporatestolerances, allowing topological differences between the target trees,and consequently making the process robust against false branches, whichare commonly observed in the segmentation (due, for example, to leaks)and skeletonization (due, for example, to naturally occurring “bumps”along the airway wall) result of in vivo data.

Reducing Computing Time

Computing the maximum clique is an NP-complete problem. For example,attempting to match two trees of 200-300 vertices each in a naïve waycould require computing time on the order of hours. Reducing the problemsize will decrease this time. Such a reduction could be achieved indifferent ways. For example, the overall problem size could be reduced,or the problem could be split into a greater number of smallersub-problems. Some embodiments of the present matching methods use acombination of both of these approaches.

In this regard, the overall problem size may be reduced by pruning outvery short (and mostly spurious) terminal branches from the input trees.Next, the input trees (e.g., two input trees) undergo a rigidregistration to bring potentially matching branch-points geometricallyclose to each other and to allow the imposition of a distancerestriction when building the association graph. With such an approach,only geometrically close branch-points are considered as possiblematches. Next, the trees may be divided into sub-trees, and only tworelatively small sub-trees may be matched at a time. A more detaileddiscussion of one manner of carrying out these general steps is nowprovided.

Pruning

Pruning short terminal branches comprises, in some embodiments of thepresent matching methods, two steps: deleting terminal branches that areshorter than a pre-defined threshold l_(th) and deleting vertices thathave only one out-edge and replacing the in-edge and out-edge with onesingle edge (stated another way, deleting vertices that have only oneremaining out-edge and connecting parent- and child-vertices directly).FIGS. 3A-3C illustrate an embodiment of these steps. FIG. 3A shows aportion of a segmented, skeletonized airway tree before pruning. FIG. 3Bshows that edge 1 is a terminal edge and has a length l<l_(th), and thesame is identified as a false branch and deleted. In FIG. 3C, the falsebranch no longer appears, and vertex A is found to have only oneout-edge. Consequently, vertex A and edges 2 and 3 are deleted and edge4 is inserted. In embodiments of the present matching methods, thepruning is applied two times in a sequence for every tree. Thus, amaximum of two airway-generations is deleted from the periphery of thetrees.

Rigid Registration

In some embodiments of the present matching methods, a rigidregistration may be performed such that the carinas of the input treesare superimposed and the angles between the corresponding main bronchiof the two trees are minimized. The carina is the first mainbranch-point of the airway tree where the trachea splits into the twomain bronchi. The carina and the left and right bronchus of each inputtree may be identified first, as follows. For every tree a depth firstsearch may be performed, which labels every edge e with the minimum andmaximum values of the x, y, and z positions of all vertices that aresituated topologically underneath e, designated as x_(min)(e),x_(max)(e), . . . , z_(max)(e). The number of vertices foundtopologically underneath e is recorded as well, designated as N(e).

The spatial extents Δ_(x)(e)=x_(max)(e)−x_(min)(e),Δ_(y)(e)=y_(max)(e)−y_(min)(e), and Δ_(z)(e)=z_(max)(e)−z_(min)(e) cannow be computed for every edge e. Next, a breadth first search may beperformed starting from the root of the tree. The carina may beidentified as the first vertex that is encountered withΔxyz=max{Δ_(x)(e), Δ_(y)(e), Δ_(z)(e)}≧50 millimeters (mm) and

$\frac{N(e)}{V} \geq 0.1$

for both of its out-edges. Similarly, the branch-points at the end ofthe two main bronchi may be found by finding the next two vertices afterthe carina that satisfy the same conditions.

The notation in the following equations refers to FIGS. 4A and 4B. Therotation angles are the x, y, and z-axes, φ_(yz avg), φ_(xz avg), andφ_(xy avg), respectively, computed with

$\begin{matrix}{{\phi_{{yz}\mspace{11mu} {avg}} = {\begin{bmatrix}{{\phi_{diff}\left( {{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TR}_{y}}^{1}}{{\overset{\rightarrow}{v}}_{{TR}_{z}}^{1}}},{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TR}_{y}}^{2}}{{\overset{\rightarrow}{v}}_{{TR}_{z}}^{2}}}} \right)} +} \\{\phi_{diff}\left( {{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TL}_{y}}^{1}}{{\overset{\rightarrow}{v}}_{{TL}_{z}}^{1}}},{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TL}_{y}}^{2}}{{\overset{\rightarrow}{v}}_{{TL}_{z}}^{2}}}} \right)}\end{bmatrix}/2}}{\phi_{{xz}\mspace{11mu} {avg}} = {\begin{bmatrix}{{\phi_{diff}\left( {{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TR}_{x}}^{1}}{{\overset{\rightarrow}{v}}_{{TR}_{z}}^{1}}},{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TR}_{x}}^{2}}{{\overset{\rightarrow}{v}}_{{TR}_{z}}^{2}}}} \right)} +} \\{\phi_{diff}\left( {{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TL}_{x}}^{1}}{{\overset{\rightarrow}{v}}_{{TL}_{z}}^{1}}},{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TL}_{x}}^{2}}{{\overset{\rightarrow}{v}}_{{TL}_{z}}^{2}}}} \right)}\end{bmatrix}/2}}{\phi_{{xy}\mspace{11mu} {avg}} = {\begin{bmatrix}{{\phi_{diff}\left( {{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TR}_{x}}^{1}}{{\overset{\rightarrow}{v}}_{{TR}_{y}}^{1}}},{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TR}_{x}}^{2}}{{\overset{\rightarrow}{v}}_{{TR}_{y}}^{2}}}} \right)} +} \\{\phi_{diff}\left( {{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TL}_{x}}^{1}}{{\overset{\rightarrow}{v}}_{{TL}_{y}}^{1}}},{\tan^{- 1}\frac{{\overset{\rightarrow}{v}}_{{TL}_{x}}^{2}}{{\overset{\rightarrow}{v}}_{{TL}_{y}}^{2}}}} \right)}\end{bmatrix}/2}}} & (3.1)\end{matrix}$

where the function φ_(diff)(φ₁, φ₂) finds the shorter of the twopossible difference angles between φ₁ and φ₂ and is defined by:

$\begin{matrix}{{\phi_{diff}\left( {\phi_{1},\phi_{2}} \right)} = \left\{ \begin{matrix}{{\phi_{1} - \phi_{2}},} & {{{if}\mspace{14mu} {{\phi_{1} - \phi_{2}}}} \leq \pi} \\{{{2\pi} = {\phi_{1} - \phi_{2}}},} & {{{{if}\mspace{14mu} \phi_{1}} - \phi_{2}} > \pi} \\{{{2\pi} = {\phi_{1} - \phi_{2}}},} & {{{{if}\mspace{14mu} \phi_{1}} - \phi_{2}} < \pi}\end{matrix} \right.} & (3.2)\end{matrix}$

The 3D transformation of the vertex coordinates can be written in thestandard matrix representation (Foley (1990)) as

[x _(T) ,y _(T) ,z _(T) ]′=T·[x,y,z]′  (3.3)

where [x, y, z] is the original vertex coordinate and [x_(T), y_(T),z_(T)] is the transformed vertex coordinate. The transformation matrix Mis defined by

$\begin{matrix}{M = {\begin{bmatrix}1 & 0 & 0 & {\overset{\rightarrow}{v}}_{T_{x}}^{1} \\0 & 1 & 0 & {\overset{\rightarrow}{v}}_{T_{y}}^{1} \\0 & 0 & 1 & {\overset{\rightarrow}{v}}_{T_{z}}^{1} \\0 & 0 & 0 & 1\end{bmatrix} \cdot R_{z} \cdot R_{x} \cdot R_{y} \cdot \begin{bmatrix}1 & 0 & 0 & {\overset{\rightarrow}{v}}_{T_{x}}^{2} \\0 & 1 & 0 & {\overset{\rightarrow}{v}}_{T_{y}}^{2} \\0 & 0 & 1 & {\overset{\rightarrow}{v}}_{T_{z}}^{2} \\0 & 0 & 0 & 1\end{bmatrix}}} & (3.4)\end{matrix}$

with

$\begin{matrix}{{R_{z} = \begin{bmatrix}{\cos \; \phi_{{xy}\mspace{11mu} {avg}}} & {{- \sin}\; \phi_{{xy}\mspace{11mu} {avg}}} & 0 & 0 \\{\sin \; \phi_{{xy}\mspace{11mu} {avg}}} & {\cos \; \phi_{{xy}\mspace{11mu} {avg}}} & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{bmatrix}}{R_{x} = \begin{bmatrix}1 & 0 & 0 & 0 \\0 & {\cos \; \phi_{{yz}\mspace{11mu} {avg}}} & {{- \sin}\; \phi_{{yz}\mspace{11mu} {avg}}} & 0 \\0 & {\sin \; \phi_{{yz}\mspace{11mu} {avg}}} & {\cos \; \phi_{{yz}\mspace{11mu} {avg}}} & 0 \\0 & 0 & 0 & 1\end{bmatrix}}{R_{y} = \begin{bmatrix}{\cos \; \phi_{{xz}\mspace{11mu} {avg}}} & 0 & {\sin \; \phi_{{xz}\mspace{11mu} {avg}}} & 0 \\0 & 1 & 0 & 0 \\{{- \sin}\; \phi_{{xz}\mspace{11mu} {avg}}} & 0 & {\cos \; \phi_{{xz}\mspace{11mu} {avg}}} & 0 \\0 & 0 & 0 & 1\end{bmatrix}}} & (3.5)\end{matrix}$

Matching

In some embodiments of the present matching methods, matching of theremaining branch-points (those other than the carina and the end-pointsof the two main bronchi) may be performed in multiple (e.g., two) stepsto help reduce computing time. For example, first only major branchpoints (which are generally branch-points that are a parent of sub-treeshaving substantial size, such as those with a spatial extent of Δxyz=30mm) are matched. This requires minimal computing time because there are,on average, only 20-30 major branch-points identified per tree. Then,starting from pairs of matched major branch-points, sub-trees arematched, only one pair of sub-trees at a time. Prior to matching twosub-trees, a simplified rigid registration may be executed bytranslating one of the two sub-trees such that the root points of thesub-trees are superimposed.

In other words, with respect to the major branch-point matching, in someembodiments of the present matching methods, major branch-points otherthan the carina and the end-points of the two main bronchi may beidentified using the same techniques discussed above with respect toidentifying the carina and the end-points of the two main bronchi duringthe rigid registration. The only difference is that the threshold forthe spatial extent is lowered to Δxyz=30 mm and no restriction for theminimal number of vertices in a sub-tree is used.

Building Association Graph

In some embodiments of the present methods, the association graph isbuilt in keeping with a goal of having its maximum clique represent asmany valid matches as possible. One could add a vertex to theassociation graph for every possible pair of vertices in the originaltrees and connect association graph vertices with edges if theyrepresent the exact same relationships between the respective verticesin the two input trees. However, considering the average number of 200to 300 vertices in the input trees, this approach could produce a hugeassociation graph, and the time to find the maximum clique will risebeyond practical values. Furthermore, such an approach does not allowfor tolerances and will not tolerate false branches well. In someembodiments of the present methods, it is beneficial to keep the size ofthe association graph small and make the overall matching processtolerant against false branches.

The segmentation and skeletonization processes may leave some false(extraneous, anatomically non-existing) branches in the tree. Treesscanned at a high lung volume are more prone to such branches than treesscanned at a low lung volume. It should not be expected that a falsebranch will simultaneously occur in both input trees (e.g., whenintra-subject trees are at issue) and at the same location.

A false branch does not change the inheritance relationship of any twovertices in the graph. Siblings stay siblings, parents stay parents oftheir respective children and vice-versa. However, the topologicaldistance between some vertices does change when a false branch isintroduced.

Vertex Relationship Array

In embodiments of the present methods, the inheritance relationship andthe topological distance between any two vertices for each input tree(e.g., two input trees) is known prior to building the associationgraph. To ascertain these properties, the two-dimensional vertexrelationship array R_(v) of dimension N×N is computed, where N standsfor the number of vertices in the associated tree. The cell R_(v)(s, t)of R_(v) contains the inheritance relationship r_(s,t) ε {PARENT, CHILD,SIBLING, N/A} and the topological distance d_(s,t) ε N between thesource vertex s and the target vertex t, with (r_(s,t)=N/A, d_(s,t)=0)

s=t. For s≠t no inheritance relationships except PARENT, CHILD, andSIBLING are possible and r_(s,t)≧1. If vertex n is not a directdescendant of vertex m and m is not a direct descendant of n thenr_(n,m)=r_(m,n)=SIBLING. No “COUSINS”, “NEPHEWS” or similar areassigned. In this regard, embodiments of the present methods and devicesmay be characterized as an automated method comprising, or a computerreadable medium comprising machine readable instructions for,respectively, (a) accessing a first representation of an airway tree ofa human subject, the first representation including vertices; (b)accessing a second representation of the airway tree of the humansubject, the second representation including vertices; and (c) buildinga tree association graph, where the building includes assigning verticesto the tree association graph such that no inheritance relationshipsother than parent, child, or sibling exist among the vertices. Usingthis inheritance model simplifies the handling of trifurcations that mayoccur in the form of two bifurcations that follow each otherimmediately, and that may occur in different order in the second tree.In either case the vertices that are topologically below thetrifurcation stay siblings and thus can still be matched. Further, thepresence of false branches does not disturb the inheritancerelationship. Consequently, it becomes possible to tolerate falsebranches to a certain degree—a benefit when matching in-vivo trees.

R_(v) may be computed with Algorithm 1 depicted in FIG. 2, which isbased on a breadth-first search. It should be noted that d_(min)=d_(max)for some embodiments of the present branch-point matching methods, andfor some embodiments of the present anatomical labeling methods it maybe d_(min)≠d_(max).

In some embodiments of the present methods, vertices are only added tothe association graph if the two corresponding vertices in the trees tobe matched are not farther than d_(max) apart. The value of d_(max) maybe set to 40 mm for matching the major branch-points, and to 15 mm formatching the sub-trees (see below).

Thus, embodiments of the present methods and devices may becharacterized as an automated method comprising, or a computer readablemedium comprising machine readable instructions for, respectively, (a)accessing a first representation of an airway tree of a human subject,the first representation including vertices; (b) accessing a secondrepresentation of the airway tree of the human subject, the secondrepresentation including vertices; and (c) building a tree associationgraph, where the building includes determining whether two possiblevertices of the tree association graph have a topological distance thatfalls within a range of acceptable topological distances.

In some embodiments of the present methods, there are guidelines for theplacement of edges in the association graph. For example, consider twovertices in the association graph, v_(assoc1) and v_(assoc2). In someembodiments, if v_(assoc1) represents a match between vertex v_(1a) intree 1 and v_(2a) in tree 2, and if v_(assoc2) represents a matchbetween vertex v_(1b) in tree 1 and v_(2b) in tree 2, then v_(assoc1)and v_(assoc2) are only connected with an edge if the following aretrue:

the inheritance relationship from v_(1a) to v_(1b) is identical to theinheritance relationship from v_(2a) to v_(2b);

the topological distance between v_(1a) and v_(1b) and the topologicaldistance between v_(2a) and v_(2b) differ by at most ∀2;

the Euclidean distance between v_(1a) and v_(1b) and the Euclideandistance between v_(2a) and v_(2b) differ by at most 20%; and

the angle between the two vectors v_(1a)-v_(1b) and v_(2a)-v_(2b) is atmost 1 radian.

These factors minimize the size of the association graph andconsequently speed up the matching process. The tolerance fortopological distances makes it possible to allow false branches to somedegree (e.g., such as the false branches not eliminated by pruning).

Finding Maximum Clique

A number of exact algorithms as well as heuristic approximationtechniques have been presented for solving the maximum clique problem(see Carraghan et al. (1990); Ostergard (2002); and Pardalos et al.(1999)). For most large problems an algorithm based on a heuristicmethods is the only reasonable choice in order to keep the computingtime within practicable limits. However, a heuristic method does notguarantee an optimal solution—something an exact algorithm can do. Inembodiments of the branch-point matching methods, the branch-pointmatching problem can be divided into a number of relatively smallsub-problems as discussed below. For that reason, an exact maximumclique algorithm is applied in certain embodiments. Algorithm 2 (basedon Carraghan et al. (1990)) in Table 1 below, which is a non-weighted(or un-weighted) maximum clique algorithm, may be used in this regard.The runtime for the matching process amounts to 1 to 3 seconds for twotrees containing 200 to 300 branch-points each (measured on a 1.2 GHzAMD Athlon™ single CPU system).

TABLE 1 Algorithm 2   Input: C: Graph with set of vertices V and set ofedges E. Output: C: Maximum unweighted clique. 1 max ← 0 2 clique ← Ø 3tmp_clique ← Ø 4 _maxClique(V, 0) 5 return clique 6 function_maxClique(U, size) 7 | if |U| = 0 then 8 | | if (size > max) then 9 | || max ← size 10 | | |_ clique ← tmp_clique 11 | |_ return 12 | while U ≠Ø do 13 | | if size + |U| ≦ max then 14 | | |_return 15 | | i ← min{j|v_(j) ∈ U} 16 | | U← U \ {v_(i)} 17 | |_ _maxClique(U ∩ neighbors(v₁),size +1) 18 | return

Implementation

In some embodiments of the present methods, the graph-related parts maybe realized with the help of the Boost Graph Library, Siek et al.(2002). The complete matching algorithm may be implemented as acommand-line application written in ANSI C++(e.g., one example ofmachine-readable instructions). Such a program can take two XML treefiles as input and can output an XML tree-matching file.

XML is a suitable platform because of the expandability of XML-basedfiles and because powerful parsers are readily available. Two differentdata files were designed:

The tree file contains all tree-related data except the originalgray-scale image and volumetric segmentation result. Information abouttopology and geometry are stored, as well as the results of allquantitative measurements. The tree file can hold trees withdisconnected parts and trees containing loops, if needed. Anatomicallabels, work status (information about what processing steps have beenperformed on the tree), and text-based notes provided by the user canalso be stored; and

The tree-matching file contains information about matchingbranch-points, segments, and lobes between two airway-trees. As with thetree file, the tree-matching file can also hold information about thework status and it can hold text notes.

C++ and Python libraries were developed for accessing these files.Providing these libraries has several advantages: data access issimplified and new applications can be added with relative ease;data-integrity is guaranteed because all read and write operations aredone through the libraries; and cross-platform operability is assured byusing standardized language features only, and by using cross-platformthird-party libraries only.

Pages 142-171 (beginning with section 8.3.4 of page 142) of Appendix 1of U.S. Provisional Patent Application Ser. No. 60/568,184, which pagesare incorporated by reference in this application, describe oneembodiment of the two file formats and the C++ and Python libraries thatcan come with them.

Modularity

The processing steps described in this disclosure may be integrated intoseparate modules. Each module may be specialized to a specific task anddata exchange may be standardized as described above. Individual modulescan be exchanged or new modules can be added as needed. Specializedtools can be built by combining different modules into a higher-levelsystem.

Modules have been developed to carry out aspects of this disclosure thatwere command-line based and parameterized through command-lineparameters. Integration into scripts may be currently be achieved via acommand-line interface (e.g., with the ‘os.system( )’ command inPython). Integration into a system based on standardized binary moduleslike the one described in pages 141-143 of Appendix 1 of U.S.Provisional Patent Application Ser. No. 60/568,184, which pages areincorporated by reference in this application, may be realized byreplacing the thin command-line interface layer with the new interfacecode.

All modules may be implemented in standard ANSI C++ for cross-platformoperability. Portable and freely available third-party libraries may beused.

Results

FIGS. 5A and 5B show a result of using an embodiment of the presentmatching methods for intra-subject matching. The run time for matchingtwo full-sized in-vivo trees (200-300 branch-points each) was about 1 to2 seconds, measured on a single AMD Athlon™ CPU running at a clock rateof 1.2 GHz.

Specifically, validation of an embodiment of the present branch-pointmatching methods was performed with in-vivo data sets. The accuracy ofthe embodiment was evaluated by comparing the automated results againstan independent standard that was based on the hand-matching resultprovided by several independent human experts.

Independent Standard

Scans from a total of 17 subjects, 10 healthy and 7 diseased, wereavailable. Each subject was scanned twice, one scan at functionalresidual capacity (FRC, 55% lung-volume), and one scan at total lungcapacity (TLC, 85% lung-volume). In each volume the airway-tree wassegmented and skeletonized by automated methods, and the graphrepresentation of the airway-trees was then available for matching. Aninteractive computer program was developed that allowed human experts toperform the tree-matching. Each pair of FRC/TLC trees was matchedindependently by three different human experts. For each tree-pair, theFRC and the TLC tree originated from the same subject. It was thereforepossible to match branch-points beyond the anatomically named points. Amatch between two branch-points was only taken as a reference if amajority of human observers (at least two out of three) agreed on it.

Automated Method

The embodiment of the present branch-point matching methods was run onall 17 tree pairs. The input trees were taken from the skeletonizationprogram “as is”; no manual pre-processing (e.g., pruning) took place.All 17 tree pairs were matched with the automated program using the samestandard parameters. No hand-adjusting of parameters took place.

For every tree pair the result of the automated matching program wasvalidated against the independent reference described in the previoussub-section. Table 2 lists the validation results for the matchingalgorithm:

TABLE 2 10 subjects 7 subjects in in vivo vivo normal diseased TotalTotal reference matches 193 211 404 Verified computer 124 130 254matches: correct 113 (91.1%) 123 (94.6%) 236 (92.9%) incorrect 11 (8.9%) 7 (5.4%) 18 (7.1%) Computer matches not 79 40 119 in referenceReference matches not 80 88 168 in computer matches

In the intra-subject branch-point matching, 92.9% of the verifiablematches agreed with the independent standard. It is interesting tocompare the 7.1% incorrect matches with the intra-observer disagreement.On average, 7.5% of the branch-point correspondences identified by thehuman experts disagreed with each other. For some tree-pairs the humanexpert disagreement was as high as 28.6%.

While the embodiments of the present matching methods (and devices) havebeen discussed in the context of the human airway tree, those ofordinary skill in the art will understand that they are also applicableto matching aspects (such as vertices) of any other artificial ornatural tree structure.

Next we turn to some embodiments of the present anatomical labelingmethods.

Assigning anatomical labels makes it possible to match airway treeacross different subjects. This allows, for example, physiologicalstudies where, across a population, the average change of airwaygeometry during the breathing cycle is of interest. Having anatomicallabels assigned to an airway tree also makes it possible to identifylobes and sub-lobes, which is crucial for surgical planning.

In some embodiments of the present methods, anatomical labels areassigned by matching the target tree against a pre-labeled tree thatrepresents a population average and contains information about thetypical geometrical and topological variability in human airway trees.Because of false branches an anatomical segment may consist of more thanone edge in the graph representation of the airway tree that is to belabeled. For that reason, labels may be assigned to the respectiveendpoint (branch-point) of the segment. An anatomical segment can haveonly one endpoint, and a branch-point can only be the endpoint of oneanatomical segment. Consequently, this approach allows for anunambiguous labeling that is independent of the number of graph-edgesthat make up one anatomical segment.

In general, the labeling of some embodiments of the present labelingmethods is performed by matching the target tree against a populationaverage. The population average may be built from a number of airwaytrees that are hand-labeled by one or more human experts. The populationaverage may incorporate the variability that can typically be observedacross subjects of the population.

Population Average

The population average is a data base that contains a reference-airwaytree with anatomical labels. It also contains the mean values forvarious geometrical measures, as well as information about the typicalgeometrical and topological variations across the population. Thefollowing sub-sections describe the details about how the populationaverage may be built according to some embodiments of the presentlabeling methods.

Notation

A vector notation may be used for computing the population average. InFIG. 6, the vector {right arrow over (s)}_(n) defines the orientationand length of segment n, and {right arrow over (s)}_(m), does the samefor segment m. Every segment has a start point, a center point, and anend point, designated by the subscripts s, c, and e, respectively. Thedifference vector {right arrow over (d)}_(nm)={right arrow over (s)}_(m)_(c) −{right arrow over (s)}_(n) _(c) is defined by the two centerpoints.

Rigid Registration

The population average will contain information about the absolutespatial orientation of segments and the spatial relation between pairsof segments. This information can only be computed relative to apredefined reference frame. The coordinates of the tree vertices aredefined in a Cartesian coordinate system. Trees used for computing thepopulation average undergo a rigid registration such that the carinalies in the origin of the coordinate system, the left main bronchus isaligned with the z-axis of the coordinate system (with the end of theleft main bronchus pointing into positive z direction), and the planedefined by the two main bronchi is congruent with the x-z plane of thecoordinate system (this also makes possible the measurement of absolutepositions of segments, though this information is not contained in theembodiment of the population average discussed because it can bereplaced with relative spatial measurements). A 3D transformation of thevertex coordinates as it is used within the embodiment of the presentbranch-point matching methods discussed above may be used.

Information Contained in Population Average

Two different classes of knowledge are contained in the populationaverage. The first class comprises pre-knowledge related to singlesegments. The second class comprises pre-knowledge about therelationship between pairs of segments.

For the single segment pre-knowledge, the following parameters arerecorded:

Segment length—for, e.g., every segment the absolute length is computed.The population average contains the mean μ_(l) and standard deviationσ_(l) of the measured values, separately listed for, e.g., everyanatomically named segment.

Spatial orientation—for, e.g., every anatomically named segment the meanvector {right arrow over (μ)}_(SO) is computed from T hand-labeled inputtrees with

$\begin{matrix}{{\overset{\rightarrow}{\mu}}_{SO} = \frac{\sum\limits_{i = 0}^{T - 1}{\overset{\rightarrow}{s}}_{n}}{T}} & (4.1)\end{matrix}$

and the unit vector of {right arrow over (μ)}_(n) is recorded in thepopulation average. Furthermore, the average variation from {right arrowover (μ)}_(n) is found by computing the standard deviation of the anglebetween {right arrow over (μ)}_(n) with

$\begin{matrix}{\sigma_{SO} = \sqrt{\frac{1}{T - 1}{\sum\limits_{i = 0}^{T - 1}\left( {\cos^{- 1}\frac{{\overset{\rightarrow}{S}}_{n} \cdot {\overset{\rightarrow}{\mu}}_{SO}}{{{\overset{\rightarrow}{S}}_{n}} \cdot {{\overset{\rightarrow}{\mu}}_{SO}}}} \right)^{2}}}} & (4.2)\end{matrix}$

FIGS. 7A and 7B show the variation of the two single segment measures,taken from 17 trees. Both the length measure and the spatial orientationmeasure have a relatively small variation and provide reliable criteriafor identifying segments.

For the inter-segment relationship the following parameters arerecorded:

Inheritance relationship—for, e.g., every pair of anatomically namedsegments the inheritance relationship is recorded as a label from theset {PARENT, CHILD, SIBLING}.

Topological distance—for, e.g., every pair of anatomically namedsegments the minimum and maximum value of the topological distance isrecorded.

Angle between segments—for, e.g., every pair {right arrow over (s)}_(n)and {right arrow over (s)}_(m) of anatomically named segments theaverage angle and its standard deviation are computed with

$\begin{matrix}{{\overset{\rightarrow}{\mu}}_{\phi} = {\frac{1}{T}{\sum\limits_{i = 0}^{T - 1}{\cos^{- 1}\frac{{\overset{\rightarrow}{S}}_{m} \cdot {\overset{\rightarrow}{S}}_{n}}{{{\overset{\rightarrow}{S}}_{m}} \cdot {{\overset{\rightarrow}{S}}_{n}}}}}}} & (4.3)\end{matrix}$

and

$\begin{matrix}{\sigma_{\phi} = \sqrt{\frac{1}{T - 1}{\sum\limits_{i = 0}^{T - 1}\left( {{\cos^{- 1}\frac{{\overset{\rightarrow}{S}}_{m} \cdot {\overset{\rightarrow}{S}}_{n}}{{{\overset{\rightarrow}{S}}_{m}} \cdot {{\overset{\rightarrow}{S}}_{n}}}} - {\overset{\rightarrow}{\mu}}_{\phi}} \right)^{2}}}} & (4.4)\end{matrix}$

respectively.

Spatial relationship between segments—for, e.g., every pair of {rightarrow over (s)}_(n) and {right arrow over (s)}_(m) of anatomically namedsegments the average spatial relationship is computed with

$\begin{matrix}{{\overset{\rightarrow}{\mu}}_{SR} = \frac{{\sum\limits_{i = 0}^{T - 1}{\overset{\rightarrow}{S}}_{m_{c}}} - {\overset{\rightarrow}{S}}_{n_{c}}}{T}} & (4.5)\end{matrix}$

and the unit vector of {right arrow over (μ)}_(SR) is recorded in thepopulation average. The average variation of the spatial relationship iscomputed as

$\begin{matrix}{\sigma_{SR} = \sqrt{\frac{1}{T - 1}{\sum\limits_{i = 0}^{T - 1}\left( {{c{os}}^{- 1}\frac{\left( {{\overset{\rightarrow}{S}}_{m_{c}} - {\overset{\rightarrow}{S}}_{n}} \right)_{c} \cdot {\overset{\rightarrow}{\mu}}_{SR}}{{{{\overset{\rightarrow}{S}}_{m_{c}} - {\overset{\rightarrow}{S}}_{n_{c}}}} \cdot {{\overset{\rightarrow}{\mu}}_{SR}}}} \right)^{2}}}} & (4.6)\end{matrix}$

FIGS. 8A-8C show the variation of various intra-segment measures. Onecan see that the angle between segments and the spatial relationshipbetween segments provide relatively good measures with only moderatevariability. On the other hand, measuring the length ratio of segmentsdoes not provide a good measure, as can be seen in FIG. 8C. Thevariability for many pairs of segments is too big.

Introducing Parallel Edges

In the embodiment of the intra-subject branch-point matching methoddescribed above, short branches (likely to be false) were pruned outfrom the target trees. While this is an appropriate and effective methodof preprocessing trees for intra-subject matching, pruning only works ina limited way when preparing a tree for anatomical labeling. FIG. 9illustrates the reason for that.

The problem is that sometimes airway trees contain relatively long“false” branches. Sometimes these branches are not actually false; theydo really exist in the target tree. They represent anatomical segmentsthat are not typical and are likely not represented in the populationaverage. The “false” branch in FIG. 9 is such an example. Pruning with afixed threshold cannot remove such branches without removing desiredbranches as well (e.g., TriLUL-LB3 in FIG. 9). For the intra-subjectmatching, the presence of such extra branches generally does not matterbecause they most likely exist in both trees. But for the labelingprocess presently described they do matter because they are notcontained in the population average and consequently they cause amismatch for the topological distances.

In some embodiments of the present labeling methods, a solution is tointroduce “parallel” edges into the target tree prior to building theassociation graph. For example, parallel edges may be introduced for alledges that might be false in the anatomical textbook-sense. For example,in FIG. 10, a parallel edge v_(p) is added to the tree if edge e₁ hastwo out-edges, among which one is a terminal edge (e₃). Additionally,the angle between e₁ and e₂ has to be close to 180 degrees. This secondcriterion my be judged by taking the ratio of distances between verticesv₁, v₂, and v₃. A parallel edge is added if |v₃−v₁|/(|e₁|+|e₂|)≧0.9 and|e₁|+|e₂|>5 voxels or |v₃−v₁|/(|e₁|+|e₂|)≧0.7 and |e₁|+|e₂|>5 voxels.

In one respect, the present methods and devices that involve theintroduction of such edges may be characterized as an automated methodcomprising, or a computer readable medium comprising machine readableinstructions for, respectively, accessing a representation of a targettree, the representation having edges and vertices; identifying a vertexin the representation that includes a first out-edge, a second out-edge,and an in-edge, where the second out-edge is a terminal edge and where afirst vertex defines one end of the first out-edge and a second vertexdefines one end of the in-edge; and introducing an edge to therepresentation that extends between the first vertex and the secondvertex. In another respect, such methods and devices may also include,or further comprise machine readable instructions for, accessing datathat includes (a) a reference tree having labels and reference treeedges, and (b) inheritance relationship data about some of the referencetree edges; and building a tree association graph, where the buildingincludes: adding a first edge to the tree association graph if (i) acorresponding edge of the representation has an inheritance relationshipthat is the same as the inheritance relationship of a correspondingreference tree edge, and (ii) the topological distance between a vertexof the first edge and a vertex of the corresponding edge of therepresentation tree is within a certain limit.

After introducing the parallel edges the mutual inheritance relationshipand topological distance for, e.g., every pair of vertices and everypair of edges may be determined. First the vertex relationship arrayR_(v) is computed using Algorithm 1 shown in FIG. 2. Based on R_(v),which contains the inheritance relationship and the minimum and maximumtopological distance between any two edges in the associated tree, theedge relationship array R_(e) is computed with Algorithm 3 shown inTable 3:

TABLE 3 Algorithm 3   Input:

: tree

: vertex relationship array Output:

: edge relationship array re-size(

, N, N) for every element in

 do | r ← “N/A” |_ d_(min) ← d_(max) ← 0 for every edge e₁t in

do | for every edge e₂t in

do | |  t₁ ← targetVertex(e₁) | |  t₂ ← targetVertex(e₂) | |  

 (e₁, e₂).r ←

 (t₁, t₂).r | | 

 (e₂, e₁).r ←

 (t₂, t₁).r | | 

 (e₁, e₂).d_(min) ←

 (t₁, t₂).d_(min) | | 

 (e₂, e₁ ).d_(min) ←

 (t₂, t₁ ).d_(min) | |  

 (e₁, e₂).d_(max) ←

 (t₁, t₂ ).d_(max) | |  

(e₂, e₁ ).d_(max) ←

 (t₂, t₁).d_(max) | |  if

 (e₁, e₂).r = “SIBLING” then | | |

 (e₁, e₂)d_(min) ←

 (e₁, e₂).d_(min)− 1 | | |

 (e₂, e₁).d_(min) ←

 (e₂, e₁ ).d_(min) − 1 | | |

 (e₁, e₂).d_(max) ←

 (e₁, e₂).d_(max) − 1 |_ |_ |_

 (e₂, e₁ ).d_(max) ←

 (e₂, e₁).d_(max) − 1Note that the relationship of every pair of edges is computed in thisexample, including parallel edges, which have a CHILD-CHILD inheritancerelationship.

Introducing parallel edges allows the embodiment of the present matchingalgorithms to chose from two possibilities—either use the original edgesand ignore the parallel edge, or use the parallel edge and ignore theoriginal edges. Two parallel edges cannot be used simultaneously in thelabeled end-result. Again referring to FIG. 10, e₁ and e_(p) have aSIBLING relationship. At the same time, e₁ is a parent of e₄ and e_(p)is a parent of e₄ as well. To make a simultaneous labeling of e₁, e_(p)and e₄ possible, all of the relationships

would have to be contained in the population average. But thisconstellation can only occur in the presence of parallel edges, andparallel edges are not allowed to be in the population average. Hence asimultaneous labeling of parallel edges is not possible.

Building Association Graph

In some embodiments of the present labeling methods, the first stepcomprises pruning the target-tree of short branches. A fixed thresholdvalue of 5 mm may be used to make sure that no essential branches arepruned away, yet obviously false branches are removed. This reduces theproblem size and shortens the overall computation time.

Next, parallel edges as described above may be added to the target-tree.

Vertices may be added to the association graph by pairing segments fromthe population average with potentially corresponding edges from thetarget tree. Edges may be added to the association graph if thecorresponding edges in the target tree have the same inheritancerelationship as the corresponding segments in the population average,and if the topological distance is within the limits given by thepopulation average.

In some embodiments of the present labeling methods, every vertex andevery edge in the association graph has a weight w_(vertex)=[0,1] andw_(edge)=[0,1], respectively, associated with it.

In some embodiments of the present labeling methods, every associationgraph vertex weight w_(vertex) is based on the single-segment measuresand is computed by

$\begin{matrix}{\omega_{vertex} = {{\exp\left( {- \frac{ - \mu_{}}{2\sigma_{l}^{2}}} \right)} \cdot {\exp\left( {- \frac{\phi_{SO}}{2\sigma_{SO}^{2}}} \right)}}} & (4.7)\end{matrix}$

with

$\begin{matrix}{\phi_{SO} = {\cos^{- 1}\frac{{\overset{\rightarrow}{S}}_{n} \cdot {\overset{\rightarrow}{\mu}}_{SO}}{{{\overset{\rightarrow}{S}}_{n}} \cdot {{\overset{\rightarrow}{\mu}}_{SO}}}}} & (4.8)\end{matrix}$

In some embodiments of the present labeling methods, every associationgraph edge weight w_(edge) is based on the inter-segment measures and iscomputed by

$\begin{matrix}{\omega_{edge} = {{\exp\left( {- \frac{\left( {\phi_{\varphi} - \mu_{\phi}} \right)^{2}}{2\sigma_{\phi}^{2}}} \right)} \cdot {\exp\left( {- \frac{\phi_{SR}^{2}}{2\sigma_{SR}^{2}}} \right)}}} & (4.9)\end{matrix}$

with

$\begin{matrix}{\phi_{\varphi} = {\cos^{- 1}\frac{{\overset{\rightarrow}{S}}_{n} \cdot {\overset{\rightarrow}{S}}_{m}}{{{\overset{\rightarrow}{S}}_{n}} \cdot {{\overset{\rightarrow}{S}}_{m}}}}} & (4.10)\end{matrix}$

and

$\begin{matrix}{\phi_{SR} = {\cos^{- 1}\frac{{\overset{\rightarrow}{S}}_{n} \cdot {\overset{\rightarrow}{\mu}}_{SR}}{{{\overset{\rightarrow}{S}}_{n}} \cdot {{\overset{\rightarrow}{\mu}}_{SR}}}}} & (4.11)\end{matrix}$

In some embodiments of the present labeling methods, only vertices andedges with a value greater than 0.1 are added to the association graph.This helps limit the size of the association graph and consequentlyreduces the computing time when searching for the maximum clique.

Finding Maximum Weighted Clique

In some embodiments of the present labeling methods, the maximum cliquein the association graph is found in a manner similar to the version ofAlgorithm 2 (Table 1 above). The only differences are that size issubstituted with weights w that were assigned to the vertices and edgesof the association graph, and lines 13 and 14 of Algorithm 2 are omittedbecause the remaining edges and vertices cannot be predicted andconsequently an early termination is not possible. The clique with thebiggest sum Ω=Σ_(i) ^(ω) _(vertex) _(i) +Σ_(j) ^(ω) _(edge) _(j) issought.

One might use the average weight of all vertices and edges as measurefor the total weight of a clique. However, this would have thedisadvantage that single edges or vertices with below-average weightwould not be included into the resulting maximum clique. Summing up theweights avoids this.

Reducing Computing Time

The problem of finding the maximum clique is NP-complete. Garey et al.(1979). One way of keeping the computing time down is to reduce theproblem size. In the case of the tree-labeling, this may be done by astepwise labeling.

In some embodiments of the present labeling methods, only a sub-tree islabeled during one labeling step. A sequence of labeling steps isillustrated in FIGS. 11A-11C. Every labeling step starts from a segmentthat already has its final label assigned to it (marked with a ‘s’ inFIGS. 11A-11C). Terminal segments (marked with a ‘t’ in FIGS. 11A-11C)are used as guides only and do not receive their final label at thispoint. The segments in-between the start segment and the terminalsegments receive their final labels.

Applying this stepwise labeling results in a significant speedup of thelabeling process. A complete tree is labeled in a few seconds, whereasan attempt of labeling an entire tree in one single step can takeseveral hours. At the same time, the accuracy of the assigned labels ispreserved.

Use Lobe-Information to Increase Quality of Labeling Result

Information about the location of lung lobes can be used to restrict thesearch space when assigning anatomical labels.

The human lung consists of five distinct lung lobes. Anatomicalpre-knowledge provides information about which anatomical airway labelsare contained in which lung lobe.

In some embodiments of the present labeling methods, given an accuratelobe segmentation, the process of assigning anatomical labels can besubdivided into five separate tasks. Each lobe is treated as a separateentity and only airways within a lobe, together with the sub-set ofanatomical labels that are known to belong to that lobe, are consideredat a time. This results in a restriction of the search space because thepossible targets where a label may be assigned to are reduced. Thisshould not only speeds up the labeling process, but more importantly itshould guarantee that every airway segment label is assigned within itscorrect lobe. As a result, the risk of mislabeling a segment should besignificantly reduced and the overall quality of the labeling resultshould increase.

Use of lobe information as just described is not implemented in theembodiment (e.g., the software) of the present labeling methodsdescribed in Appendix 1 of U.S. Provisional Patent Application Ser. No.60/568,184.

Implementation

In some embodiments of the present labeling methods, the graph-relatedparts may be realized with the help of the Boost Graph Library, Siek etal. (2002). The complete labeling algorithm may be implemented as acommand-line application written in ANSI C++ (e.g., one example ofmachine-readable instructions). Such a program can take one XML treefile as input and can write the labels back into the same file. Adescription of relevant XML files appears above. Labeling a full-sizedin-vivo tree containing 200-300 vertices requires a computing time ofabout 5 seconds, measured on a single AMD Athlon™ CPU running at a clockrate of 1.2 GHz.

In some embodiments of the present labeling methods, a separatecommand-line application computes the population average. A text filedefines the trees (contained in XML tree files) that can be used for thedefinition of the population average. The results are saved into anASCII file (about 80 kBytes in size) and will be loaded by the labelingprogram during initialization. The ASCII format allows platformindependent usage of the population average while being one of the moststraightforward formats to define and implement. Computing thepopulation average may take up to a half a minute per input tree(depending on the tree size), but may be done only once.

Results

FIG. 12 shows a detail view of a typical result of automated anatomicallabeling of an in-vivo tree according to one embodiment of the presentlabeling methods. Note that correct labels are assigned despite thepresence of false branches.

Validation of the embodiment of our branch-point labeling programpresented here (minus the use of lobe-information described above) wasperformed with in-vivo data sets. The accuracy of the program wasevaluated by comparing the automated results against an independentstandard provided by a human expert.

Independent Standard

Scans from a total of 17 subjects, 10 healthy and 7 diseased, wereavailable. For each subject, a total lung capacity scan (TLC, 85%lung-volume) was taken as the basis for generating the independentstandard. In each volume the airway-tree was segmented and skeletonizedby automated methods, and the graph representation of the airway-treeswas then available to the labeling and matching processes.

An application was written that allowed human experts to perform thetree-labeling by hand. Using this tool, all 17 TLC trees werehand-labeled by a human expert. These hand-labeled trees were then usedas gold-standard for the validation of the automated method

Automated Method

The leave-one-out (jackknife) method was used for building thepopulation average and for testing the automated labeling. All trees butone were used for building the population average, and the automatedprogram (the program may, as a whole, also be characterized as analgorithm) was then run on the one tree left out. This procedure wasrepeated for all 17 TLC trees.

Tables 4 and 5 list the validation results for the labeling algorithm.In the first case (Table 4) the segment labels were verified foraccuracy. In the second case (Table 5) the branch-point labels werevalidated.

TABLE 4 Computer label Reference label Correct Incorrect not found inthe not found among Subject label label reference set computer labelsNormal 1 31 1 0 6 (brp001tlc) Normal 2 28 0 4 7 (brp002tlc) Normal 3 320 0 7 (brp003tlc) Normal 4 29 0 2 10 (brp004tlc) Normal 5 18 2 7 11(brp005tlc) Normal 6 28 0 5 5 (brp006tlc) Normal 7 30 1 1 10 (brp007tlc)Normal 8 32 0 0 9 (brp008tlc) Normal 9 25 1 7 4 (brp009tlc) Normal 10 222 8 8 (brp010tlc) Diseased 1 23 2 2 6 (h005tlc) Diseased 2 27 1 1 7(h006tlc) Diseased 3 27 0 0 10 (h008tlc) Diseased 4 33 0 0 5 (h009tlc)Diseased 5 24 2 1 10 (h012tlc) Diseased 6 24 2 4 9 (h014tlc) Diseased 728 0 2 7 (h016tlc) Total 461 14 (97.1%) (2.9%)

TABLE 5 Computer label Reference label Correct Incorrect not found inthe not found among Subject label label reference set computer labelsNormal1 29 3 0 6 (brp001tlc) Normal2 27 1 4 7 (brp002tlc) Normal3 29 3 07 (brp003tlc) Normal4 26 3 2 10 (brp004tlc) Normal5 15 5 7 11(brp005tlc) Normal6 22 6 5 5 (brp006tlc) Normal7 24 7 1 10 (brp007tlc)Normal8 29 3 0 9 (brp008tlc) Normal9 24 2 7 4 (brp009tlc) Normal10 21 38 8 (brp010tlc) Diseased 18 7 2 6 (h005tlc) Diseased 24 4 1 7 (h006tlc)Diseased 23 4 0 10 (h008tlc) Diseased 28 5 0 5 (h009tlc) Diseased 21 5 110 (h012tlc) Diseased 19 7 4 9 (h014tlc) Diseased 25 3 2 7 (h016tlc)Total 404 71 (85.1%) (14.9%)85.1% of the assigned labels are correct. The number of correctlyassigned segment labels is higher because it is difficult for theautomated algorithm to label terminal branch-points. A terminalbranch-point is the highest-generation branch-point with anatomicallabel within a branch. Because no more labels are to be found beyond aterminal branch-point, the labeling algorithm does not have guidesegments when assigning the label. This makes it difficult to find theexact branch-point. FIG. 13 illustrates this issue.

It was surprising and unexpected that the automated labeling identifieda number of human errors. After reviewing the cases of disagreementbetween the independent standard and the automated result, we were ableto identify 5 cases where a label was misplaced by the human expert.

While the embodiments of the present labeling methods (and devices) havebeen discussed in the context of the human airway tree, those ofordinary skill in the art will understand that they are also applicableto labeling aspects of any other artificial or natural tree structurefor which a naming convention exists or is desired.

The present methods and devices are not intended to be limited to theparticular forms disclosed. Rather, they are to cover all modifications,equivalents, and alternatives falling within the scope of the claims.The appended claims are not to be interpreted as includingmeans-plus-function limitations, unless such a limitation is explicitlyrecited in a given claim using the phrase(s) “means for” and/or “stepfor,” respectively.

REFERENCES

Each of the following references is specifically incorporated byreference in each location it is cited above:

-   Ballard and Brown, In: Computer Vision. Prentice Hall PTR, 1982.-   Boyden, In: Segmental anatomy of the lungs, McGraw-Hill, 1955.-   Carraghan and Pardalos, Operations Research Letters, 9:375-382, 11    1990.-   Foley, Computer Graphics: Principles and Practice, Addison-Wesley,    1990.-   Garey and Johnson, In Computers and intractability, a guide to the    theory of NP-completeness, W. H. Freeman and Company, 1979.-   Kitaoka et al., In: Automated Nomenclature Labeling of the Bronchial    Tree in 3D-CT Lung Images, Tokyo, Japan, 1-11, 2002.-   Mori et al., IEEE Transactions on Medical Imaging, 19:103-114, 2000.-   Ostergard, Discrete Applied Mathematics, 195-205, 2002-   Palágyi et al., In: Proc. 18th Int. Conf. Information Processing in    Medical Imaging, IPMI2003, Ambleside, UK, Springer, 222-233, 2003.-   Pardalos et al., In: High Performance Algorithms and Software in    Nonlinear Optimization, HPSNO97 conference Ischia, Italy, June,    297-300, 1999.-   Park, In: Registration of linear structures in 3-D medical images,    Osaka Univers., Japan, Dept. of Informatics and Mathematical    Science, 2002.-   Pelillo et al., IEEE Transactions on Pattern Analysis and Machine    Intelligence, 21:1105-1120, 1999.-   Pisupati et al., In: Mathematical Methods in Biomedical Image    Analysis, Proc. of the Workshop, 160-169, 1996a.-   Pisupati et al., In: Geometric Tree Matching with applications to 3D    Lung Structures, Proc. of the 12^(th) Annual Symposium on    Computational Geometry, Philadelphia, Pa., 419-20, 1996b.-   Siek et al., In: The Boost Graph Library, Addison-Wesley, 2002.-   Tschirren et al., In: Branchpoint Labeling and Matching in Human    Airway Trees, SPIE Medical Imaging, San Diego, Calif., 2003.-   Tschirren et al., In: Radiological Society of North America, RSNA,    88th Scientific Assembly and Annual Meeting, Chicago, Ill., 2002a.-   Tschirren et al., In: Segmentation, Skeletonization, and Branchpoint    Matching—A Fully Automated Quantitative Evaluation of Human    Intrathoracic Airway Trees, MICCAI, Tokyo, Japan, 2002b.

1. A computer readable medium comprising machine readable instructionsfor: (a) accessing a first representation of an airway tree of a humansubject; (b) accessing a second representation of the airway tree of thehuman subject; and (c) automatically pruning at least a portion of thefirst representation.
 2. The computer readable medium of claim 1,further comprising machine readable instructions for: (d) automaticallypruning at least a portion of the second representation.
 3. The computerreadable medium of claim 2, the first and second representations of theairway tree each including vertices and one or more terminal branches,and (c) including: deleting from the first representation of the airwaytree any terminal branch having a length shorter than a pre-definedthreshold length.
 4. The computer readable medium of claim 3, where (d)includes: deleting from the second representation of the airway tree anyterminal branch having a length shorter than a pre-defined thresholdlength.
 5. The computer readable medium of claim 4, where (c) alsoincludes: deleting from the first representation of the airway tree anyvertex having an in-edge and only one out-edge; and replacing thein-edge and out-edge of each deleted vertex with a single edge.
 6. Thecomputer readable medium of claim 6, where (d) also includes: deletingfrom the second representation of the airway tree any vertex having anin-edge and only one out-edge; and replacing the in-edge and out-edge ofeach deleted vertex from the second representation with a single edge.7. The computer readable medium of claim 1, the first and secondrepresentations of the airway tree each including a carina and two mainbronchi, each main bronchi including an end point, the computer readablemedium further comprising machine readable instructions for: identifyingthe carina and the end points of the two main bronchi in the first andsecond representations of the airway tree.
 8. The computer readablemedium of claim 7, further comprising machine readable instructions for:matching the identified carina and end points of the two main bronchifrom the first representation with the identified carina and end pointsof the two main bronchi from the second representation.
 9. The computerreadable medium of claim 8, further comprising machine readableinstructions for: identifying some remaining major branch-points in thefirst and second representations of the airway tree, the remaining majorbranch-points comprising branch-points that correspond to an entry pointof either a lobe or a sub-lobe.
 10. The computer readable medium ofclaim 9, further comprising machine readable instructions for: matchingthe identified remaining major branch-points from the firstrepresentation with the identified remaining major branch-points fromthe second representation.
 11. The computer readable medium of claim 10,further comprising machine readable instructions for: matching asub-tree from the first representation with a sub-tree from the secondrepresentation, where the sub-tree from a given representation is belowthe identified remaining major branch-points from the givenrepresentation.
 12. The computer readable medium of claim 11, where eachsub-tree includes a root point, and prior to matching the sub-tree fromthe first representation with the sub-tree from the secondrepresentation, one of the two sub-trees is translated such that theroot points of the two sub-trees are superimposed.
 13. The computerreadable medium of claim 1, further comprising machine readableinstructions for: building a tree association graph.
 14. The computerreadable medium of claim 13, where the tree association graph has amaximum clique, and the computer readable medium also comprises machinereadable instructions for: identifying the maximum clique.
 15. Thecomputer readable medium of claim 13, where the tree association graphincludes vertices, and the computer readable medium also comprisesmachine readable instructions for: assigning the vertices to the treeassociation graph such that no inheritance relationships other thanparent, child, or sibling exist among the vertices.
 16. The computerreadable medium of claim 1, where the first representation of the airwaytree is based on a volumetric scan of the human subject taken usingcomputed tomography.
 17. The computer readable medium of claim 1, wherethe first representation of the airway tree is based on a volumetricscan of the human subject taken using magnetic resonance.
 18. Thecomputer readable medium of claim 1, where the airway tree is diseased.19. A computer readable medium comprising machine readable instructionsfor: (a) accessing a first representation of an airway tree of a humansubject, the first representation including vertices; (b) accessing asecond representation of the airway tree of the human subject, thesecond representation including vertices; and (c) building a treeassociation graph, where the building includes determining whether twopossible vertices of the tree association graph have a topologicaldistance that falls within a range of acceptable topological distances.20. A computer readable medium comprising machine readable instructionsfor: (a) accessing a first representation of an airway tree of a humansubject, the first representation including vertices; (b) accessing asecond representation of the airway tree of the human subject, thesecond representation including vertices; and (c) building a treeassociation graph, where the building includes assigning vertices to thetree association graph such that no inheritance relationships other thanparent, child, or sibling exist among the vertices.